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SUMMARY 


We are developing an efficient, flexible, / 'Mree- 
dimensional, hydrodynamic, computer code for a spherical cap 
geometry. The code will be used to simulate NASA's Atmos- 
pheric General Circulation Experiment ( AGCE ) . The AGCE is a 
spherical, baroclinic experiment which will model the large- 
scale dynamics of our atmosphere; it has been proposed to 
NASA for future Spacelab flights. In the AGCE a radial 
dielectric body force will simulate gravity, with hot fluid 
tending to move outwards. In order that this force be 
dominant, the AGCE must be operated in a low gravity environ- 
ment such as Spacelab. 

The full potential of the AGCE will only be real- 
ized by working in conjunction with an accurate computer 
model. Proposed experimental parameter settinqs will be 
checked first using model runs. Then actual experimental 
results will be compared with the model predictions. This 
interaction between experiment and theory will be very 
valuable in determining the nature of the AGCE flows and 
hence their relationship to analytical theories and actual 
atmospheric dynamics. 

Although much is understood about the averaged 
steady-state behavior of our atmosphere, most of the time- 
dependent activity is still unexplained. We expect, from 
previous experimental work with cyclindrical baroclinic 
laboratory models, that the AGCE will also exhibit many 
periodic and irregular flows. It is the major thrust of the 
AGCE program to analyze these strongly nonlinear, unsteady 
flows . 


The computer model is based on the complete set of 
fully nonlinear, incompressible equations. Previous computer 
codes are inadequate for the accurate simulations required. 
It is known, from theory, axisymmetric computations and the 
cylindrical experiments, that the AGCE type flows have thin 
thermal layers and plumes. This necessitates nonuniform 
finite-difference meshes, and makes the hydrostatic assump- 
tion a questionable approximation. Further, the wide range 
of flow types and parameters, and the resulting large number 
of cases, imply a requirement for efficient state-of-the-art 
coding, with implicit methods, able to utilise the new fifth- 
generation vector processors. 

The model is being made flexible so that it can be 
used for design studies for other possible future spherical 
Spacelab experiments dealing with the oceans or other plane- 
tary atmospheres, and for studies in cylindrical qeometry. 
The model will be properly validated using available accurate 
flow measurements. 

The development, testing and application of this 
code is a substantial task. We originally proposed to 
develop this three-dimensional code in 1981, in a three year 
program. Funding constraints have so far allowed us to 
perform only about 40% of the work, in the first two years. 
We estimate that it will take a further three years to 
complete. During the third year the code and diagnostics 
will be completed and tested, and cylindrical validations 
will be performed. During the fourth year spherical valida- 
tions will be performed, with timing and convergence tests. 
During the fifth year we will model a variety of proposed 
AGCE configurations and tests. 


Section 1 
BACKGROUND 


1.1 THE CYLINDRICAL ANNULUS FLOWS 

One of the major unsolved problems of meteoroloqy 
is the persistence and irregularity of our atmosphere. 
Although the mean steady general circulation is quite well 
understood, the time-dependent behavior is not. The larqe- 
scale atmospheric flow sometimes forms a relatively stable 
pattern that persists for many days or weeks, while at other 
times it is much more transitory. 

A fundamental attack on these problems has been 
made through the study of well-defined laboratory model 
flows. In particular, since the ISSO's a substantial amount 
of experimental work has been devoted to the study of baro- 
clinic flows in a rotating and differentially heated cy- 
lindrical annulus of liquid (References 1’ through 5). 
Systematic scaling of the governing equations for both our 
atmosphere and the cylindrical annulus reveals that the 
annulus is a model of the large-scale (synoptic) atmosphere 
flow (Reference 6). As the parameters in the experiment are 
varied, axisymmetric flows, steady wave flows, periodically 
fluctuating wave flows (vacillations) and irreqular flows are 
observed. Figure 1 is a schematic diaqram of the cylindrical 
annulus, Figure 2 shows a typical wave flow and Figure 3 is a 
flow regime diagram. 

Analytical studies have contributed to an under- 
standing of these baroclinic flows, but owing to the mathe- 
matical difficulties associated with cylindrical geometry, 
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p - General point having cylindrical polar 
coordinates (r,0,z) in frame rotating 
with the apparatus 

«fi =(0,0,£i)- rotation vector 

g =(0,0,“9)“ acceleration of gravity 

b,a t d - fluid occupies region 

a< r<b, 0<z<d [l+n 2 (r 2 -l/2(b 2 4a 2 ))/2gd]“d 

T(r,0,z,t) - temperature at general point P and time t 

Ta ,Tb - T(a,0,z,t) ; T(b,0 t z,t) respectively 

Figure 1 Schematic Diagram Illustrating Rotating 
Fluid Annulus 
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Figure 2 Photograph of a Typical Annulus Wave Flow 
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Figure 3 Experimental results on the transitions between 
the flow regimes and wave numbers for a fixed 
depth of water in the annulus and for fixed radii 
of the cylinders, plotted in the ^4 versus ^ 2 
diagram. Experimental d_etails: a = ~3.48cm, b = 

6.02cm, d = 10.00cm, T = 20.0°C, 1.01.x 

10" 2 cm 2 sec" :L , P = 0.998gm cm -3 , ^5 = 7.19. 
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non-separability and nonlinearity, these studies have dealt 
mainly with simplified models and linearized stability 
problems (References 7, 8 and 9). Numerical studies have 
dealt with the basic axisymmetric states and with a steady 
wave flow (References 10, 11 and 12). Much has been learned 
from the annulus work about baroclinic instability and baro- 
clinic waves, and insight into the behavior of larqe-scale 
atmospheric dynamics has been achieved (References 13 and 
14 ). 


1.2 THE AGCE 

* 

The major inadequacy of the cylindrical annulus 
configuration as a model of the atmosphere is the absence of 
spherical curvature. A true spherical model requires a 
radial gravity field. 


In 1978, a proposal for a spherical baroclinic 
experiment known as the Atmospheric General Circulation 
Experiment (’AGCE) was submitted to the Global Weather 
Program, NASA/OSSA. In the apparatus proposed, the liquid is 
held between two concentric spheres, with a larqe AC voltaqe 
applied across them. The resultinq temperature-dependent 
dielectric body force corresponds exactly to the inward 
gravity acceleration, 


_ 2 eyV 2 r r t r b } 2 
E pa l r t - r 1 


( 1 ) 


Here e and p are the dielectric constant and the density, y 
and a are the corresponding thermal coefficients, rb and 
rt are the inner and outer radii, V is the voltage and r is 
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the radius. In order for gg to be dominant the AGCE 
apparatus will have to be flown in an orbitinq vehicle such 
as Spacelab. Figure 4 is a schematic drawing of a proposed 
apparatus. This figure also includes some preliminary esti- 
mates for the dimensions and imposed conditions. Tempera- 
tures are maintained on the spherical boundaries such that 
the liquid experiences a stable radial gradient and a lati- 
tudinal gradient. The spheres are rigidly corotated. 
Scientific and engineering design studies for the AGCE are 
under way. 

Substantial differences between the AGCE flows and 
the cylindrical annulus flows are anticipated. It is the 
major objective of the AGCE program to study the non-linear 
unsteady flows occurring in the spherical model. 

Another spherical geophysical fluid dynamics 
experiment which also exploits the radial dielectric body 
force was proposed earlier than the AGCE and will be flown on 
Spacelab 3. This experiment is known as the Geophysical 
Fluid Flow Cell (GFFC) and is concerned with convection due 
to radially unstable temperature gradients. The scientific 
teams for the AGCE and the GFFC are listed in Table 1. 

1.3 The AGCE DESIGN PROGRAM 

The major objective of the AGCE is to study non- 
linear, baroclinic wave flows in a spherical geometry. Thus 
the experiment must allow a parameter range in which axisym- 
metric flows are strongly baroclinically unstable. A number 
of preliminary analytical design study calculations for the 
AGCE have been performed and the results indicate that to 
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Figure 4 Schematic Drawing of the Proposed AGCE Apparatus 











Table 1 

SCIENTIFIC TEAMS FOR SPACELAB EXPERIMENTS 


Geophysical Fluid Flow Cell ( GFFC ) 


Principal Investigator: John Hart (Univ. of Colorado) 

Co-investigators: Juri Toomre (Univ. of Colorado) 

Peter Gilman ( NCAR/HAO) 

George Fichtl ( NASA/MS FC ) 
William Fowlis (NASA/MSFC) 

Fred Leslie (NASA/MSIC) 


Atmospheric General Circulation Experiment ( AGCE) 


Principal Investigator: William Fowlis (NASA/MSFC) 

Co-investigators: George Fichtl (NASA/MSFC) 

John Geisler (Univ. of Utah) 
Robert Gall (Univ. of Arizona) 
Basil Antar (Univ. of Tennessee 
Space Institute) 

Glyn Roberts (SAI) 

Eric Pitcher (Univ. of Miami) 
Timothy Miller (NASA/MSFC) 



obtain strong instability the apparatus will have to be care- 
fully designed (References 15 through 23). However, because 
of the simplifications which have to be made to perform 
analytical calculations, these results must be considered 
qualitative in nature and cannot be used to prepare specifi- 
cations for the apparatus. For quantitative theoretical 
work, numerical models are required. 

1.4 THE AXISYMMETRIC AND STABILITY CODES 

A numerical modeling design program for the AGCE is 
under way (References 31 through 34). The approach taken was 
to develop two numerical models. The first model obtains 
steady, axisymmetric basic states for the spherical qeometry 
and is based on the complete, fully nonlinear set of equa- 
tions. The second model examines the stability of these 
basic states. It is based on the complete linearized set of 
pertubation equations. We are running these two-dimensional 
models for systematically varied values of the AGCE dimen- 
sions, fluid properties and imposed conditions. This work is 
producing theoretical regime diagrams comparable to Fiqure 3, 
which can be used to prepare specifications for the AGCE 
apparatus such that strong baroclinic instability will be 
achieved. 


We validated these models by applying them to a 
cylindrical geometry. Experimental and numerical data are 
available for the basic states of the cylindrical annulus 
flows, and experimental data are available for the regime 
diagrams (References 3, 4, 10, 11, 12, 24 and 25). Excellent 
agreement has been obtained. We (have also compared the pre- 
dictions of the stability code with the experimental regime 
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diagrams. Again, the agreement was excellent. Finally, we 
validated the spherical basic state model using laser-Doppler 
measurements of spin-up in spherical geometry (cf., Reference 


Section 2 

THREE-DIMENSIONAL NUMERICAL MODEL FOR THE AGCE 

2.1 MODEL OBJECTIVES AND JUSTIFICATION 

Although a substantial amount of work with two- 
dimensional numerical models is under way (see Section 1.4), 
the full scientific potential of the AGCE will not be 
realized unless a three-dimensional, fully nonlinear, numeri- 
cal model based on the complete set of equations is also 
available. We are developing such a model. 

The model has two primary objectives. First, it 
will demonstrate the possible types of nonlinear baroclinic 
wave flows that may occur in AGCE for different parameter 
ranges. This will aid in the selection of the sequence of 
experimental parameters for the AGCE on Spacelab. Secondly, 
the numerical results will extend the limited measurements 
that can be made with the AGCE apparatus. Owing to the 
relatively small and confined volume of fluid and to the 
presence of the large AC electric field, it will not be 
possible to measure the fluid behavior thoroughly. Techni- 
ques for flow and temperature measurement for the AGCE have 
been proposed (References 15 and 27). If the measurements 
agree with the model predictions, this will tend to confirm 
both the model predictions of the whole flow and temperature 
field, and also the model predictions for other cases. 



2.2 


SPECIFICATION OF THE PROBLEM 


2.2.1 Model Equations 

We are using the Boussinesg convection equations in 
the spherical polar coordinates (0,X,r), with corresponding 
velocity components u, v and w. Significant components of 
our system of equations are summarized in Table 2. The 
divergence of the stress and momentum tensor has many extra 
terms associated with the curvilinear coordinate system? 
these are handled correctly. 

The gravity potential $ has up to three components, 
as shown in Tables 2 and 3. The radial component corresponds 
to an inward gravity field Gt.(r/r|-)P; the radial powers 
of particular interest are -2 for atmospheric gravity and -5 
for the gravity field which models the effect of the di- 
electric body force, cf., Equation (1). 

Axial (terrestrial) gravity is included for simu- 
lating laboratory experiments in a cylindrical annulus and in 
prototype AGCE configurations. Centrifugal gravity is 
included since theory and annulus experiments both demon- 
strate that it is very significant for the larqer rotation 
rates. It is absent in the atmosphere, since the vertical 
coordinate is based on equipotent ial surfaces with the 
centrifugal potential included. 

2.2.2 Computational Domain and Special Cases 

The axisymmetric computational domain in spherical 
polar coordinates is defined by the inequalities 


EQUATIONS FOR THE MODEL 


Component 

Description 

Continuity 

Incompressible 

Formulation 

Primitive variable, not hydrostatic 

Primary variables 

T, u, v, w, p 

Thermal Diffusivity 

k constant 

Kinematic viscosity 

v constant 

Buoyancy acceleration 

bV$ 

Specific volume increase 

b = aT + a 2 T 2 

Heat equation 

Optional internal sources or sinks 

Coriolis force 

Coordinates rotate at R 

Gravity 

Radial, axial and centrifugal 


Table 3 

GRAVITY POTENTIAL COMPONENTS 


+ 


Component 

Formulation 

Radial 

Terrestrial 

Centrifugal 

9 t r P+1 /(P+1)r t P 

gz 

-R 2 p 2 /2 
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e* < e < e r » 

r b < r < r t» 
and shown in Figure 5. 



Figure 5 General Form of the Axisymmetric Domain 


Special cases of particular interest are listed 
in Table 4. The extensive set of cylindrical annulus 
experimental results will provide an important means of code 
validation. The annulus 

P* < P < P r , 

0 < z < h . 
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Table 4 

SPECIAL CASES FOR THE COMPUTATIONAL DOMAIN 


Case 

Values 

Spherical layer 
Hemisphere Cap 
Cylindrical Annulus 

e * - 0 ' e r - 

9* ■ 0 , e r = /2 

~ P A /R ' e r ~ p r^ R 

r b » R r fc - R + h 


can be approximated by using a large value of R. An inter- 
mediate value of 9 r is required for the spherical layer of 
Figure 4. 


2,2.3 Options for the Boundary Conditions 

There are two basic options for the temperature 
boundary conditions on each of the four boundary segments, 

T - imposed distribution, or 

V ■ 0 

The second option is appropriate for a thermally insulatinq 
boundary. Naturally, at least one boundary must have an 
imposed temperature distribution. 

The four imposed boundary temperature distributions 
are determined by interpolation from imposed corner temper- 
atures, T T br , T fcJl , and T tr . Possible values are 15°C, 


2-5 



25°C, 25 °C and 35 # C. The top and bottom temperatures are 
linear in an imposed power of 0, in cos 0/a or in In 

(tan 1 0/2), at the user's option. The side temperatures are 
linear in an imposed power of r. Other options could easily 
be added later as the need develops. Examples are to impose 
the normal heat flux, or to impose the heat flux as a 
function of the local temperature. 

There are four options for the flow boundary condi- 
tions on each of the four boundary segments. The first is 
the no-slip condition, with the solid boundary rotatinq, 
relative to the rotating coordinate system, at R s - Q. The 
second is the free-slip condition appropriate to a free 

surface (and to the equator if it is a symmetry boundary). 
The third is the axis condition. These three options are 

summarized in Table 5. The fourth simulates an unresolved 

turbulent 'sublayer by making the stress linear in the slip 
velocity; this option is intended for atmosphere simulations. 

Table 5 

VELOCITY BOUNDARY CONDITION OPTIONS 


Boundary Type 

Boundary Conditions 

Location 

No-slip 

c 

u 

£ 

u 

o 

■ 


v = (fi - Q)r sin 0 

All 

Free-slip 

u = a Q v = a e w=o 

V e r 


w = a r v = a r u = o 

r b ' r t 

Axis 

o 

ii 

CD 

CO 

11 

> 

11 

a 

O 

11 

°< 

CD 



0 = % 



r 
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2.2.4 


The X Boundaries 


We assume periodicity in the X direction, with 
period 2%/m^. This is naturally always valid with m^ * 1. 
By using larger values we can obtain higher resolution at 
much lower cost. 

2.3 SPATIAL REPRESENTATION 

We use second order finite differences on nonuni- 
form meshes to represent the 0 and r derivatives, as in the 
present two-dimensional codes (References 31 through 34). We 
presently use finite differences for the X derivatives: we 
have not used the attractive fast Fourier transform algorithm 
because of its poor vectorizat ion properties on the Cyber 
205. However, we have written the code in such a way that 
the Fast Fourier Transform option could be added at a later 
stage. 


2.3.1 The Nonuniform Meshes 


The mesh spacings for the 0 and r meshes are pro- 
portional to the product of the distances from points just 
outside the domain. This gives small spacings in the 
boundary layers without wasting mesh points in the interiors. 
Our X mesh is uniform. 

2.3.2 Placement of the Variables 


We use a staggered computational mesh. The primary 
variables (T and p) are collocated, and are not defined on 
the boundary but on points half a mesh interval inside and 
outside. The velocity components u, v and w are defined at 
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the natural locations for representing the corresponding 
components of the pressure gradient as differences of adja- 
cent pressure values. This staggered placement of the vari- 
ables requires averaging in the representation of the 
Coriolis and buoyancy forces, and slightly increases the 
total amount of spatial averaging required in the model. 
However, it gives improved accuracy in the representation of 
the pressure gradient and the velocity divergence, and makes 
the pressure equation (a Poisson representation) much more 
simple than with the non-staggered mesh in our original 
proposal . 


2*3.3 Conservation Properties 

Our spatial representation conserves representa- 
tions of total heat and angular momentum, and conserves 
representations of the T 2 and kinetic energy integrals under 
the action of the advection, Coriolis force, and pressure 
gradient terms. These are all desirable properties for any 
spatial representation of these equations, since they help 
ensure stability and aid in the diagnostics and interpreta- 
tion. 


2.4 TEMPORAL REPRESENTATION 

2.4.1 Time-Stepoing and Iterating to a Steady Solution 

The code has two objectives and two corresponding 
modes of operation. One is to determine the detailed time 
evolution by ordinary time stepping. The other is to iterate 
to a steady solution. In the second case, we use an itera- 
tion method developed for our two-dimensional AGCE studies 
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(References 31-34) in which the time steps can be different 
for each mesh point, and different for temperature and 
velocity. Thus large steps can be taken in the interior, 
while smaller steps are taken near the boundary where the 
mesh spacing is smaller. 

2.4.2 Alternating Direction Implicit Algorithm 


We use an ADI treatment for advection, diffusion 
and Coriolis effects. This allows large time steps, which is 
a very important consideration in a three-dimensional code. 
These algorithms were also used in the two-dimensional AGCE 
studies and in the earlier three-dimensional ocean fore- 
casting code (References 31-36). 

2.4.3 Internal Waves 


We use a leapfrog representation of internal 
gravity waves, so that the temperature and velocity fields 
are updated alternately. This representation is subject to 
the stability criterion N6t < 2, which can be restrictive due 
to large N values in thin thermal layers. 

We optionally evade this stability criterion by 
modifying the short wavelength components of the change in 
temperature, in the layers where N 2 is large. Again, this 

method has already been used successfully in ocean prediction 
(Reference 36). We used a simpler algorithm in the two 
dimensional work, where the detailed time dependence was not 
of primary concern. 
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2.5 THE PRESSURE ALGORITHM 

2.5.1 Treatment of the Pressure Gradient 

The pressure field is determined by the continuity 
equation and the boundary condition of zero normal flow. We 
obtain a first approximation to the new velocity field by 
performing an implicit time step (or iteration) using an 
initial pressure field obtained by extrapolation from the 
total pressure field at the two previous time steps. 

We then correct the velocity field, so that its 
divergence is zero and there is no normal flow, by adding the 
gradient of a correction pressure field. This correction 
pressure is determined by solving a representation of 
Poisson's equation. 

This pressure algorithm is used in our two- 
dimensional AGCE models (References 31 through 34). It was 
used previously in the codes described in References 45 and 
46. 

2.5.2 Solution of Poisson Equation 

The pressure equation is solved by using the fast 
Fourier transform in the \ direction and by using ADI itera- 
tion with Jordan's optimum parameters in the transverse 
plane. Convergence is rapid, and this method does not 
increase computation time unduly. Two Fourier transforms are 
required at each time step. If significantly better Poisson 
methods for our application become available, we will use 
them. 
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2.6 

INPUT, OUTPUT, 

DIAGNOSTICS AND GRAPHICS 




The problem, 

numerical 

method , 

and ou 

tput 

are 

defined 

and controlled 

using a 

data page 

in a 

standard 

format . 

This allows 

us to use 

default 

values 

for 

most 


parameters, and to calculate a sequence of cases in one run. 

Resulting fields are output to disc files at 
regular intervals. Thi.® allows us to use separate graphics 
programs without repeating the calculations, and to split 
long computations into multiple segments. 

Significant diagnostic quantities (maxima, inte- 
grals, etc.) are printed out in tabular form every few time 
steps or iterations. This will allow us to monitor conver- 
gence and detect when and where problems occur. 

Printer graphics are used for testing and develop-- 
ment, to provide rapid display of results. Other more 
sophisticated forms of graphical display will be developed 
and used for reporting purposes. 

2.7 CONVERGENCE AND ACCURACY 

♦ 

We will analyze the accuracy of our steady and 
time-dependent solutions by repeating some of our computa- 
tions, using Smaller time steps and more mesh intervals. 
This will provide an accurate measure of convergence. 
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2.8 


CODING CONSIDERATIONS AND VECTOR HARDWARE 


We intend that this code should be as portable as 
possible, consistent with reasonable speed. We anticipate 
use on several different computers, by various NASA employees 
and contractors. We will therefore write most of the code in 
FORTRAN 66. 


The speed constraint is associated with the vector 
hardware on modern computers. To make full use of this 
hardware on the Cyber 205 requires a certain amount of non- 
standard FORTRAN. 

We hope eventually to automate the process of 
transferring the code from computer to computer. Editor 
programs will • translate the non-standard features from 
dialect to dialect. We have usd such editors previously to 
transfer a code between IBM and CDC non-vector computers and 
the vectorizing Texas Instruments ASC {Reference 44). 


Section 3 

VALIDATION OP THE CODE 


3.1 SPHERICAL GEOMETRY 

Although it will not be possible to validate the 
three-dimensional code in a spherical geometry until the AGCE 
apparatus is flown in Spacelab, some limited validation is 
possible. Our present plans are described in the next two 
paragraphs . 


First, consider the results obtained in the numeri- 
cal design program for the theoretical AGCE regime diagram 
(see Section 1.3). If we use the basic states as initial 
conditions and then insert small-amplitude wave perturba- 
tions, we should obtain agreement with the previously derived 
theoretical regime diagram. 

Secondly, the code will be able to treat special 
spherical cases for which laboratory experiments can be 
performed and measurements obtained. These include the AGCE 
apparatus, together with various preliminary test versions, 
because our code includes the option of terrestrial (axial) 
gravity. They may also include unstratified spherical flows, 
with differential rotation producing non-axisymmetric motion. 

3.2 CYLINDRICAL GEOMETRY 

Since the three-dimensional code will be designed 
with the flexibility for easy conversion to cylindrical 
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geometry, it can be checked against known annulus wave flows 
(References 2S-30). 

The annulus experiments have demonstrated that as 
the parameters (temperature difference and rotation rate) are 
varied, a steady three-dimensional solution with (for 
example) four waves can break down, and then evolve in a 
nonlinear way to a five-wave steady solution. It should be 
possible to simulate this instability and thus to validate 
the code. 


The annulus experiments have also demonstrated 
periodic flows (vacillations). First, there is amplitude 
vacillation, in which a four-wave solution fluctuates between 
almost axisymmetric surface flow and a strongly asymmetric 
flow. Secondly, there is wave number vacillation, with the 
flow cycling in a complex manner between a mainly four-wave 
structure and a mainly five-wave structure. We would hope to 
be able to simulate both types of vacillation. Related 
vacillations are of course a major feature of atmospheric 
flows . 


Section 4 

OTHER APPLICATIONS 


In addition to its value for the AGCE proqram, the 
i hree-d imensional code will have several other applications. 

In its cylindrical form, we expect to qain new 
understanding of significant remaining theoretical problems 
for the annulus time-dependent flows. 

The code could also be modified to study ideas for 
other future, AGCE-type, Spacelab model experiments, simu- 
lating perhaps ocean or mantle flows, or stellar or planetary 
convection . 

The code could be modified to simulate aspects of 
convection flows in a planetary or stellar interior, within 
the limits of the Boussinesq approximation. 

The Boussinesq approximation could be modified to 
allow a mean density dependent on the total gravitational 
potential (the anelastic approximation). This would broaden 
the range of as trophysical problems to which it could be 
applied . 
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